(*  :Title:	Inverse z-Transforms  *)

(*  :Authors:	Brian Evans, James McClellan  *)

(*
    :Summary:	One-dimensional and separable multidimensional
		inverse z-transforms
 *)

(*  :Context:	SignalProcessing`Digital`InvZTransform`  *)

(*  :PackageVersion:  2.7	*)

(*
    :Copyright:	Copyright 1989-1991 by Brian L. Evans
		Georgia Tech Research Corporation

	Permission to use, copy, modify, and distribute this software
	and its documentation for any purpose and without fee is
	hereby granted, provided that the above copyright notice
	appear in all copies and that both that copyright notice and
	this permission notice appear in supporting documentation,
	and that the name of the Georgia Tech Research Corporation,
	Georgia Tech, or Georgia Institute of Technology not be used
	in advertising or publicity pertaining to distribution of the
	software without specific, written prior permission.  Georgia
	Tech makes no representations about the suitability of this
	software for any purpose.  It is provided "as is" without
	express or implied warranty.
 *)

(*  :History:	*)

(*  :Keywords:	z-transform, region of convergence  *)

(*  :Source:	*)

(* 
    :Warning:	Keep conditional clauses on the same line; separating them
		  by RETURNS will confuse Mathematica (e.g., the package
		  context is never reset to Global); you can place a
		  newline after an operator (like &&)

		The n argument of myinvz must always be a symbol NOT
		  an expression (see postinvzrules and local shift operator)
 *)

(*  :Mathematica Version:  1.2 or 2.0  *)

(*  :Limitation:  *)

(*
    :Discussion:  1-D Local state rules base has 62 rules:
		    I.   multidimensional hooks     2 rules
		    I.   rational transforms       26 rules
		    II.  non-rational transforms    9 rules
		    III. transform properties      12 rules
		    IV.  transform DSP structures   4 rules
		    V.   transform strategies       9 rules

	Note that  the  partial  fractions  decomposition  strategy is
	implemented  as  two  rules --   one in the rational transform
	pairs section and one in the properties section.

	At each step in the inverse rules base, the current expression
	has a local state associated with it.   This state consists of
	a list of nine boolean values.   Each boolean value is associ-
	ated with a strategy.   If  an  element  is  True,  then  that
	strategy has not been tried yet;  if False, then that strategy
	has already been tried, and it will not be tried again.  Thus,
	local state is used  to  prevent infinite loops which would be
	caused by the repetitive application of certain strategy rules.
	See the section  S T A T E  D E F I N I T I O N  below and see
	section V of the rules.

	Unlike the  other  symbolic  transforms, these rules are main-
	tained in a list called InvZTransformRules.   This  was neces-
	sary because Mathematica  reordered these rules in undesirable
	ways if they were coded  as a set of transformation functions.
	Another benefit of  the list form is that it allowed more con-
	trol over how rules are applied.

	In general,  this rule base works on expressions with negative
	powers of z.  Great  pains  are  taken  to  convert  any input
	expressions  to  this form accordingly.  This conversion works
	very well for rational z-transform functions.

	The  InvZTransform  function  first  calls the interface rules
	(InvZTransformInterfaceRules), which handle default arguments,
	error messages, and multidimensional transforms.   These rules
	simply call  MyInvZTransform  once per dimension of the trans-
	form. This inverse z-transform is biased toward causal systems.
	However,  by specifying the proper region of convergence (ROC),
	anti-causal sequences can be found (see InvZTransform::usage).

	The driver for the one-dimensional  rule base is the six argu-
	ment version of MyInvZTransform.  MyInvZTransform[X, z, n, rm,
	rp, s, op] either returns the  inverse transform as a function
	of n or an approximation to  the actual inverse transform.  If
	the rule base can do neither,  which  should  not  happen, the
	interface function cleanup will report any errors.   Arguments
	of MyInvZTransform:

	X  z-transform function		rm  Rminus of ROC
	z  z-transform variable		rp  Rplus of ROC
	n  "time"-domain variable	st  local state semaphores
					op  options

	The data tag myinvz is  manipulated by the InvZTransformRules:
	myinvz[X, z, n, rm, rp, st, op, zlist, rmlist, rplist].    The
	last three arguments are only used  by the special multidimen-
	sional rules.
 *)

(*  :Functions:	InvZTransform  *)



(*  B E G I N     P A C K A G E  *)

BeginPackage[ "SignalProcessing`Digital`InvZTransform`",
	      "SignalProcessing`Digital`ZSupport`",
	      "SignalProcessing`Support`TransSupport`",
	      "SignalProcessing`Support`ROC`",
	      "SignalProcessing`Support`SigProc`",
	      "SignalProcessing`Support`FilterSupport`",
	      "SignalProcessing`Support`SupCode`" ]


If [ TrueQ[ $VersionNumber >= 2.0 ],
     Off[ General::spell ];
     Off[ General::spell1 ] ];


(*  U S A G E     I N F O R M A T I O N  *)

InvZTransform::usage =
	"InvZTransform[e,z] or InvZTransform[e,z,n] gives the inverse \
	z-transform of e. \
	Here, e can be a function of z, a list, or a z-transform object \
	(see ZTransData). \
	As a list, e should contain three elements:  F[z], rminus, and rplus, \
	such that the region of convergence is rminus < |z| < rplus. \
	If e is a z-transform object, then the inverse z-transform is simply \
	InvZTransform[e] because e contains all necessary information. \
	This routine returns signal processing expressions. \
	Convert them to mathematical formulas by using TheFunction."

(*  E N D     U S A G E     I N F O R M A T I O N  *)


Begin["`Private`"]


(*  U S E R     I N T E R F A C E  *)

(*  InvZ operator  *)
Unprotect[InvZ]
InvZ/: TheFunction[ InvZ[z_, n_][f_] ] := InvZTransform[f, z, n]
Protect[InvZ]

(*  InvZTransform  *)
InvZTransform/: Options[InvZTransform] :=
	{ Apart -> Rational, Dialogue -> True, Simplify -> True,
	  Terms -> 10, TransformLookup -> {} }

InvZTransform[x_] :=
	invzdriver[ Options[InvZTransform], x ]
InvZTransform[x_, z_] :=
	invzdriver[ Options[InvZTransform], x, z ]
InvZTransform[x_, z_, n_] :=
	invzdriver[ Options[InvZTransform], x, z, n ]
InvZTransform[x_, z_, n_, op__] :=
	invzdriver[ ToList[op] ~Join~ Options[InvZTransform], x, z, n ]

(*  Global variables  *)
completeOptions = {}
dialogueFlag = False
dialogueAllFlag = False
indentationString = ""
simplifyFlag = False

nList = {}
zList = {}
multidFlag = False
rmList = {}
rpList = {}

(*  validvarQ  *)
validvarQ[ x_Symbol ] := True
validvarQ[ x_List ] := Apply[And, Map[VariableQ, x]]
validvarQ[ x_ ] := False

(*  Driver for inverse z-transform interface  *)
invzdriver[op_, x_, z_, n_] :=
	Message[ Transform::badvar, "z", InvZTransform, z ] /;
	! validvarQ[z]

invzdriver[op_, x_, z_, n_] :=
	Message[ Transform::badvar, "discrete-time", InvZTransform, n ] /;
	! validvarQ[n]

invzdriver[op_, args__] :=
	Block [ {},
		dialogue = Replace[Dialogue, op];

                (* Set all global variables right here *)
		completeOptions = op;
		dialogueFlag = TrueQ[dialogue] || SameQ[dialogue, All];
		dialogueAllFlag = SameQ[dialogue, All];
		simplifyFlag = TrueQ[Replace[Simplify, op]];
		nList = zList = {};
		multidFlag = False;

                (* Do the transform  *)
		cleanup [ invztransform[op, args], InformUserQ[op] ] ]

invztransform[op_, args__] :=
	Replace [ zinverse[op, args], InvZTransformInterfaceRules ]

InvZTransformInterfaceRules = {
	zinverse[op_, e_] :>
		invztransform[ op, e, ZVariables[e] ] /; ZTransformQ[e],
	zinverse[op_, e_] :>
		Message[ InvZTransform::novariables, "z", GetVariables[e] ],

	zinverse[op_, e_, z_] :>
		invztransform[op, e, z, DummyVariables[Length[z], Global`n]],

	zinverse[op_, e_List, z_, n_] :>
		invztransform[op, MakeZObject[e, z], z, n],
	zinverse[op_, e_, z_Symbol, n_] :>
		MyInvZTransform[e, z, n, op],
	zinverse[op_, e_, z_List, n_] :>
		multiDInvZTransform[e, z, n, op]
}

cleanup[ Null, flag_ ] := Null

cleanup[ trans_, flag_ ] :=
	Block [	{},
		If [ flag, explain[trans]; Scan[explain, trans, Infinity] ];
		trans ]

explain[ myinvz[f_, z_, rest__] ] :=
	Message[ Transform::incomplete, "inverse z-transform", f, z ]


(*  M E S S A G E S  *)

InvZTransform::badterms = "Non-integer terms for series expansion: ``"
InvZTransform::badROC = "Improper region of convergence in ``."


(*  T R A N S F O R M A T I O N     S T A T E  *)

	(*  Listed in order of appearance in the rule base  *)

polyfactorfield = 1		(* factor denominator of polynomials	     *)
expcheckfield = 2		(* check to see if r^n form		     *)
normdenomfield = 3		(* normalize denominator		     *)
partialfractionsfield = 4	(* partial fractions on poly w Apart         *)
partialfractionsfield2 = 5	(* partial fractions on poly w MyApart	     *)
logformfield = 6		(* factoring inside of a logarithmic form    *)
expformfield = 7		(* factoring inside of an exponential form   *)
expandallfield = 8		(* apply ExpandAll[] to current expression   *)
anticausalfield = 9		(* try anti-causal expansion		     *)
cepstrumfield = 10		(* apply rule for complex cepstrum to log    *)
seriesexpansionfield = 11	(* power series expansion of cur. expression *)

statevariables = 11		(* number of state variables *)

(*
    return the initial transformation state, which is a list of whether or not
    partial fractions, the first denominator normalization, and the second
    denominator initialization has been performed yet on that expression.
 *)

initInvZstate[] := Table[True, {statevariables}]
nullInvZstate[] := Table[False, {statevariables}]
anticausalInvZstate[] :=
	Block [	{state},
		state = SetStateField [	nullInvZstate[], logformfield, True ];
		SetStateField [ state, expformfield, True ] ]

resetstate[s_] := SetStateField[s, expcheckfield, True]


(*  S U P P O R T I N G     R O U T I N E S  *)

(*    Since the inverse z-transform rule base is biased toward causal	*)
(*  sequences, the leftOrRightSided rules adjust the right-sided	*)
(*  inverse returned by InvZTransform.					*)
(*    One of the strategies handles the general case.			*)

leftsided[p_, n_] :=
	( p /. { Step[b_. n + c_.] :> - Step[- b n - b + c] } )
rightsided[p_, n_] := p

leftOrRightSided[a_, p_, n_, rm_, rp_] := leftsided[p, n]  /; SameQ[Abs[a], rp]
leftOrRightSided[a_, p_, n_, rm_, rp_] := leftsided[p, n]  /; SameQ[Abs[-a], rp]
leftOrRightSided[a_, p_, n_, rm_, rp_] := leftsided[p, n]  /; N[Abs[a] > rp]
leftOrRightSided[a_, p_, n_, rm_, rp_] := rightsided[p, n] /; SameQ[Abs[a], rm]
leftOrRightSided[a_, p_, n_, rm_, rp_] := rightsided[p, n] /; SameQ[Abs[-a], rm]
leftOrRightSided[a_, p_, n_, rm_, rp_] := rightsided[p, n] /; N[Abs[a] < rm]
leftOrRightSided[a_, p_, n_, rm_, rp_] := rightsided[p, n] /; InfinityQ[rp]
leftOrRightSided[a_, p_, n_, rm_, rp_] := rightsided[p, n]

(* downsampledRewrite *)
downsampledRewrite[ f_, k_, M_, z_ ] :=
	Block [ {},
		downsampledExpRewrite[f, k, M, z] /.
		  ( z^b_. :> takeRealPower[z^b, M] ) ]

(* downsampledExpRewrite  *)
downsampledExpRewrite[ f_, k_, M_, z_ ] :=
	Block [ {},
		f /. ( z^a_. Exp[c_ k / M] :> z^a /;
			 FreeQ[a, z] && (c == -a M I Pi 2) ) ]

(* downsampledQ  *)
downsampledQ[ f_, k_, M_, z_ ] :=
	(! FreeQ[f, z]) && (! FreeQ[f, M]) &&
	  (! SameQ[ f, downsampledExpRewrite[f, k, M, z] ])

(*  fixUp -- strip redundant options and remove TransformLookup option  *)
fixUp[ op_ ] :=
	{ Dialogue -> Replace[Dialogue, op],
	  Simplify -> Replace[Simplify, op],
	  Terms -> Replace[Terms, op] }

(*  mDresamplingCheck  *)
mDresamplingCheck[Upsample, x_, z_] :=
	Block [ {cond, upmatrixtrans, zvars},
		cond = ! MyFreeQ[ GetAllFactors[x, z], zList ];
		If [ cond,
		     Needs[ "SignalProcessing`Digital`MDInvZTransform`" ];
		     cond = upsampleSystemQ[x, zList] ];
		cond ]

(*  multByExponential  *)
multByExponentialQ[f_, z_] :=
	Block [	{scale},
		scale = ScalingFactor[f, z];
		( ! SameQ[scale, 1] ) && ( ! SameQ[scale, 0] ) ]

(*  multiDInvZTransform  *)
multiDInvZTransform[e_, z_, n_, op_] :=
	Block [	{invztrans, oldflag, oldnList, oldrmlist, oldrplist, oldzList},
		{oldzList, oldnList, oldFlag, oldrmlist, oldrplist} =
			{zList, nList, multidFlag, rmList, rpList};
		zList = z;
		nList = n;
		multidFlag = True;
		If [ ZTransformQ[e],
		     {rmList, rpList} = { GetRMinus[e], GetRPlus[e] },
		     {rmList, rpList} = { {}, {} } ];
		invztrans = MultiDInvTransform [ e, z, n, op, ZTransformQ,
						 MyInvZTransform, MakeZObject,
						 Global`n ];
		{zList, nList, multidFlag, rmList, rpList} =
			{oldzList, oldnList, oldFlag, oldrmlist, oldrplist};
		invztrans ]

(*  normalizeDenominator  *)
normalizeDenominator[ratpoly_, z_] :=
	Block [	{adjdenom, adjnumer, const, denom, maxexp, numer, term},
		numer = Numerator[ratpoly];
		denom = Expand[Denominator[ratpoly]];
		maxexp = Exponent[denom, z];
		term = z^(-maxexp);
		adjdenom = Distribute[term denom];
		const = ConstantTerm[adjdenom, z];
		adjdenom = Distribute[ (1/const) adjdenom ];
		adjnumer = If [ MixedPolynomialQ[numer, z],
				Distribute[ term Expand[numer] ],
				term numer ];
		1/const adjnumer / adjdenom ]

(*  partialFractionsDialogue  *)
partialFractionsDialogue[p_, partfrac_, options_, apart_] :=
	Block [ {},
		If [ InformUserQ[options],
		     Print["( after the term"];
		     Print["  ", p];
		     Print["  has been broken up into its partial " ];
		     Print["  fractions representation using ", apart];
		     Print["  ", partfrac, " )" ] ] ]

(*  partialFractionsQ  *)
partialFractionsQ[p_, z_] :=
	MixedPolynomialQ[ Expand[Denominator[p]], z ] &&
	  ! FreeQ[ Denominator[p], z ]

(*  rootROC  *)
SetAttributes[rootROC, {Listable}]
rootROC[x_, y_] := takeRealPower[x, y]

(*  seriesDialogue  *)
seriesDialogue[x_, expansion_, terms_, options_] :=
	Block [	{},
		If [ InformUserQ[options],
		     Print["( after breaking up the expression"];
		     Print["  ", x];
		     Print["  into its series representation"];
		     Print["  ", expansion, " )"] ] ]

(*  seriesExpansion  *)
seriesExpansion[ True, x_, z_, power_ ] := Series[x, {z, 0, power}]

seriesExpansion[ False, x_, z_, power_ ] :=
	Block [	{negx, expansion},
		negx = x /. z -> z^-1;
		expansion = Series[negx, {z, 0, power}];
		expansion /. z -> z^-1 ]

(*  seriesStrategy  *)
seriesStrategy[ x_, z_, n_, rm_, rp_, s_, op_, rest___ ] :=
	Block [ {expansion, exponents, posexpandflag, state, terms},
		terms = Replace[Terms, op];
		state = SetStateField[s, seriesexpansionfield, False];

		(*  If the Terms option is disabled, then take no action.  *)
		If [ ! terms,
		     Return[ myinvz[x, z, n, rm, rp, state, op, rest] ] ];

		If [ ! IntegerQ[terms],
		     Message[ InvZTransform::badterms, terms ];
		     myinvz[x, z, n, rm, rp, state, op, rest] ];

		(*  See if we should expand in pos. or neg. powers	*)
		exponents = GetAllExponents[x, z];
		posexpandflag = If [ Apply[And, Thread[exponents >= 0]],
				     True, False, False ];

		(*  First expansion  *)
		expansion = seriesExpansion[posexpandflag, x, z, terms - 1];

		(*  Second expansion if first failed  *)
		If [ ! SameQ[Head[expansion], SeriesData],
		     expansion =
		       seriesExpansion[Not[posexpandflag], x, z, terms - 1] ];

		(*  If this fails, call MyInvZTransform again		*)
		(*  but with series expansion disabled			*)
		If [ ! SameQ[Head[expansion], SeriesData],
		     myinvz[ x, z, n, rm, rp, state, op, rest ],
		     seriesDialogue[x, expansion, terms, op];
		     state = nullInvZstate[];
		     Map [ myinvz[#1, z, n, rm, rp, state, op, rest]&, 
			   Normal[expansion] ] ] ]

(*  takeRealPower --  assume that both arguments are real-valued  *)
SetAttributes[takeRealPower, {Listable}]
takeRealPower[ 0, r_ ] := 0
takeRealPower[ Infinity, r_ ] := Infinity
takeRealPower[ x_^b_, r_ ] := x^(b r) /; IntegerQ[b r]
takeRealPower[ x_, r_ ] := x^r

(*  upsampledZtransform  *)
upsampledZtransform[ f_, z_, n_, rm_, rp_, s_, op_, rest___] :=
	Block [	{newf, newrm, newrp, upindex},
		upindex = UpsampleFactor[f, z];
		If [ IntegerQ[upindex],
		     upindex = Abs[upindex],
		     If [ SameQ[Replace[Dialogue, op], All],
			  upindex = Replace[upindex, -k_ :> k];
		          Print[ "assuming ", upindex, " is an integer" ] ] ];
		newf = f /. ( z^a_. :> z^(a/upindex) );
		newrm = rootROC[rm, upindex];	 (* adjust ROC *)
		newrp = rootROC[rp, upindex];
		Upsample[upindex, n]
			[myinvz[newf, z, n, newrm, newrp, s, op, rest]] ]


(*  B E G I N     R U L E     B A S E  *)

ROCrewriterules = {
	b_?Positive / Abs[a_] :> Abs[b / a],
	b_?Negative / Abs[a_] :> - Abs[- b / a]
}

(*  MyInvZByPass--- called by inverse DTFT and interfaces to MyInvZTransform  *)
MyInvZByPass[ e_, z_, n_, zlist_, nlist_, op_ ] :=
	MyInvZTransform[ e, z, n, op ]

(*  MyInvZTransform  *)
MyInvZTransform[ e_, z_, n_, op_ ] :=
	Block [ {rminus, rplus, valid},
		rminus = SPSimplify[ GetRMinus[e] ];
		rplus = SPSimplify[ GetRPlus[e] ] /. ROCrewriterules;
		valid = Apply[And, Thread[ToList[rminus] < ToList[rplus]]];
		If [ TrueQ[! valid],
		     Message[ InvZTransform::badROC, e ],
		     MyInvZTransform[ TheFunction[e], z, n, rminus,
		     		      rplus, initInvZstate[], op ] ] ] /;
	ZTransformQ[e]

MyInvZTransform[ e_, z_, n_, op_ ] :=
	MyInvZTransform[ e, z, n, 0, Infinity, initInvZstate[], op ]


(*   Driver for one-dimensional rule base			*)
(*   First, convert all z /(z - a) forms to 1/(1 - a z^-1)	*) 
(*	since the rule base favors terms with z^-1 terms	*)
(*   Loop until the expression to be inverted does not change.	*)
MyInvZTransform[ e_, z_, n_, rm_, rp_, st_, options_ ] :=
	Block [ {laste, newe, newinvzrules, nvars, op, trace},

                (* generate the z-transform rules *)
		newinvzrules = TransformFixUp[ z, n, options, myinvz, False,
					       InvZTransform, Null, Null ];
		InvZTransformRules = Join[ newinvzrules,
					   OriginalInvZTransformRules ];

		(* strip redundant options and set dialogue level *)
		op = fixUp[options];
		trace = SameQ[ Replace[Dialogue, op], All ];

		(* take the 1-D inverse transform *)
		newe = myinvz[ MapAll[convert[#, z]&, e],
			       z, n, rm, rp, st, op, zList, rmList, rpList ];
		While [ ! SameQ[ laste, newe ],
			If [ trace, Print[ newe ]; Print[ "which becomes" ] ];
			laste = newe;
			newe = laste /. ( myinvz[a__] :>
				 Replace[myinvz[a], InvZTransformRules] ) ];

		(* fix up the result *)
		newe = laste /. postinvzrules;
		While [ ! SameQ[ laste, newe ],
			If [ trace, Print[ newe ]; Print[ "which becomes" ] ];
			laste = newe;
			newe = laste /. postinvzrules ];
		If [ trace, Print[ newe ] ];

		(* simplify the resulting expression if requested *)
		If [ Replace[Simplify, op],
		     nvars = If [ EmptyQ[nList], {n}, nList ];
		     newe = SPSimplify[ newe,
					Dialogue -> trace,
					Variables -> nvars ] ];

		newe ]

(*  normalize as much of the expression as possible first  *)
convert[ z_^k_ ( z_ + a_ )^j_, z_ ] := ( 1 + a z^-1 )^j  /; j == -k
convert[ z_ / ( z_ + a_ ), z_ ] := 1 / ( 1 + a z^-1 )
convert[ x_, z_ ] := x

postinvzrules = {
	shift[a_. t_, n_, n0_] :>
		If [ Count[{t}, f_[p___][b___], Infinity] > 0,
		     a Shift[n0, n][t],			(* operators present *)
		     a (t /. n -> n - n0) ] /;		(* no operators *)
		FreeQ[a, n],

	(*  These two rules are needed to support a series	*)
	(*  expansion which returns a constant term plus	*)
	(*  terms of the form  approx z^r  .			*)
	myinvz[ b_, z_, n_, rest___ ] :>
		b Impulse[n] /;
		FreeQ[b, z],		   (* Converges all ROC *)

	myinvz[ a_. z_^r_., z_, n_, rest___ ] :>
		a Impulse[n + r] /;
		FreeQ[{a,r}, z],	   (* Converges all ROC *)

	(*  Attempt a series expansion about z = 0		*)
	(*  This is the strategy when all else has failed.	*)
	(*  The user can disable this rule (Terms -> False).	*)
	(*  Please keep this as the last rule.			*)
	myinvz[x_, z_, n_, rm_, rp_, s_, op_, rest___] :>
		seriesStrategy[x, z, n, rm, rp, s, op, rest ] /;
		GetStateField[s, seriesexpansionfield]
}

Format[ myinvz[ x_, z_, n_, rest___ ] ] := 
	SequenceForm[ ColumnForm[{"Z" Superscript[-1],
				  "  " ~StringJoin~ ToString[z]}],
		      { x } ]
Format[ shift[t_, n_, n0_] ] := SequenceForm[ {t}, Subscript[n -> n + n0] ]


(*  Most of the rest of the file is the list of inverse z-transform rules  *)

InvZTransformRules = {}

OriginalInvZTransformRules = {


  (*  M U L T I D I M E N S I O N A L     H O O K S  *)

  myinvz[ f_, z_, n_, rm_, rp_, s_, op_, rest___ ] :>
	upsampleSetUp[ f, z, n, rm, rp, s, op, rest, nList ] /;
	multidFlag && mDresamplingCheck[Upsample, f, z],

  myinvz[ newROC[ f_, args___ ], z_, n_, rest___ ] :>
	newROC[ myinvz[ f, z, n, rest ], args ],




  (*  I.  R A T I O N A L     T R A N S F O R M S		*)


  (*      A.  Simple transform pairs				*)
  (*	      1.  Constant (handles b=0 case as well)		*)
  myinvz[ b_, z_, n_, rest___ ] :>
	b Impulse[n] /;
	FreeQ[b, z],					(* Converges all ROC *)

  (*	          Redundant rule --  placed here for speed	*)
  myinvz[ a_. z_^r_., z_, n_, rest___ ] :>
	a Impulse[n + r] /;
	FreeQ[{a,r}, z],				(* Converges all ROC *)

  (*	      2.  Discrete-time sine functions			*)
  myinvz[ 1 / (1 - 2 Cos[w0_] z_^-1 + z_^-2), z_, n_, rm_, rp_, rest___ ] :>
	leftOrRightSided[1, Sin[w0 + w0 n] Step[n] / Sin[w0], n, rm, rp] /;
	FreeQ[w0, z],

  myinvz[ 1 / (1 + c_. z_^-1 + z_^-2), z_, n_, rm_, rp_, rest___ ] :>
	Block [	{sinw0, w0},
		w0 = ArcCos[- c / 2];
		sinw0 = Sqrt[1 - c^2/4];
		leftOrRightSided[ 1, Sin[w0 n + w0] Step[n] / sinw0,
				  n, rm, rp ] ] /;
	N[-2 < c < 2],

  (*	      3.  Discrete-time cosine functions		*)

  myinvz[ (1 - Cos[w0_] z_^-1) / (1 - 2 Cos[w0_] z_^-1 + z_^-2), z_, n_,
	  rm_, rp_, rest___ ] :>
	leftOrRightSided[1, Cos[w0 n] Step[n], n, rm, rp ] /;
	FreeQ[w0, z],

  (*	      4.  Discrete-time hyperbolic functions		*)

  myinvz[ (Sinh[b_] + Sinh[c_] z_^-1) / (1 - 2 Cosh[w_] z_^-1 + z_^-2), z_, n_,
	  rm_, rp_, rest___ ] :>
	leftOrRightSided[1, Sinh[b + w n] Step[n], n, rm, rp] /;
	FreeQ[{w,b,c},z] && ( c == w - b ),

  myinvz[ (Cosh[b_] + Cosh[c_] z_^-1) / (1 - 2 Cosh[w_] z_^-1 + z_^-2), z_, n_,
	  rm_, rp_, rest___ ] :>
	leftOrRightSided[1, Cosh[b + w n] Step[n], n, rm, rp] /;
	FreeQ[{w,b,c},z] && ((c == w - b) || (c == b - w)),


  (*      B.  First order sections					*)
  (*	      most covered by multiplication-by-exponential property	*)
  (*	      1/(1+az^-1) is more general than 1/(1-az^-1)		*)

  (*	      1.  Discrete-time pulse has no poles			*)
  myinvz[ ( 1 + b_. z_^L_ ) / ( 1 + a_. z_^-1 ), z_, n_, rm_, rp_, rest___ ] :>
	Block [ {},
		Assuming[ L < 0, dialogueAllFlag ];
		If [ dialogueAllFlag, Print[ "where ", L, " is an integer" ] ];
		(-a)^n Pulse[-L, n] ] /;
	FreeQ[{a,b,L}, z] && Implies[ NumberQ[L], IntegerQ[L] && L < 0 ] &&
	  ( Equal[a, b, 1] || SameQ[-b, (-a)^(-L)] ),

  (*	      2.  General exponential times a step (right-sided pair)	*)
  myinvz[ 1 / ( 1 + a_. z_^-1 ), z_, n_, rm_, rp_, rest___ ] :>
	leftOrRightSided[a, (-a)^n Step[n], n, rm, rp] /;
	FreeQ[a,z],

  (*	      3.  General first-order section				*)
  myinvz[ a_. / (c_ + d_. z_^-1), z_, n_, rm_, rp_, rest___ ] :>
	leftOrRightSided[d/c, a/c (-d/c)^n Step[n], n, rm, rp] /;
	FreeQ[{a,c,d}, z],

  myinvz[ (a_. + b_. z_^-1) / (c_ + d_. z_^-1), z_, n_, rm_, rp_, rest___ ] :>
	leftOrRightSided[d/c,
			 a/c (-d/c)^n Step[n] + b/c (-d/c)^(n-1) Step[n-1],
			 n, rm, rp] /;
	FreeQ[{a,b,c,d}, z],


  (*      C.  Second order sections					*)

  myinvz[ z_^-1 / ( 1 + b_. z_^-1 + a_. z_^-2 ), z_, n_, rm_, rp_, rest___ ] :>
	Block [	{costheta, fun, rho, theta},
		rho = Sqrt[a];
		costheta = -b / ( 2 rho );
		theta = ArcCos[costheta];
		fun = ( 1 / ( rho Sin[theta] ) ) rho^n Sin[n theta] Step[n];
		leftOrRightSided[Sqrt[a], fun, n, rm, rp] ] /;
	N[(a > 0) && (4 a > b^2)],

  myinvz[ (1 + b_. z_^-1) / (1 + c_. z_^-1 + a_. z_^-2), z_, n_,
	  rm_, rp_, rest___ ] :>
	Block [	{costheta, fun, rho, theta},
		rho = Sqrt[a];
		costheta = - b / rho;
		theta = ArcCos[costheta];
		fun = rho^n Cos[n theta] Step[n];
		leftOrRightSided[Sqrt[a], fun, n, rm, rp] ] /;
	N[ (c == 2 b) && (0 < a < b^2) ],

  myinvz[ 1 / ( 1 + b_. z_^-1 + a_. z_^-2 ), z_, n_, rm_, rp_, rest___ ] :>
	Block [ {costheta, rho, sintheta, theta},
		rho = Sqrt[a];
		costheta = - b / ( 2 rho );
		theta = ArcCos[costheta];
		sintheta = Sin[theta];
		fun = rho^n ( Cos[n theta] + Cot[theta] Sin[n theta] ) Step[n];
		leftOrRightSided[Sqrt[a], fun, n, rm, rp] ] /;
	N[ (a > 0) && (4 a > b^2) ],

  myinvz[ 1 / ( 1 + a_. z^-2 ), z_, n_, rm_, rp_, s_, op_, rest___ ] :>
	Block [ {},
		Assuming[ Negative[a], dialogueAllFlag ];
		leftOrRightSided[ Sqrt[-a],
				  ( Sqrt[-a]^n + (- Sqrt[-a])^n ) Step[n] / 2,
				  n, rm, rp ] ] /;
	FreeQ[a,z] && Implies[ NumberQ[N[a]], N[a < 0] ],

  myinvz[ 1 / ( 1 + a_. z_^-2 ), z_, n_, rm_, rp_, s_, op_, rest___ ] :>
	Block [ {},
		Assuming[ Positive[a], dialogueAllFlag ];
		leftOrRightSided[ Sqrt[a],
				  Sqrt[a]^n Cos[n Pi / 2] Step[n],
				  n, rm, rp ] ] /;
	FreeQ[a,z] && Implies[ NumberQ[N[a]], N[a > 0] ],

  myinvz[ 1 / ( 1 + a_. z_^-1 )^2, z_, n_, rm_, rp_, rest___ ] :>
	leftOrRightSided[a, (n + 1) (-a)^n Step[n], n, rm, rp] /;
	FreeQ[a,z],


  (*      D.  Higher order sections					*)
  (*	      Note that k can be real-valued.				*)

  (*	      1.  General denominator section to an integer power	*)
  myinvz[ c_. (1 + b_. z_^-1)^k_, z_, n_, rm_, rp_, rest___ ] :>
	leftOrRightSided[b, c Binomial[n-k-1,-k-1] (-b)^n Step[n], n, rm, rp] /;
	FreeQ[{b,c,k}, z] && ( k < -1 ),

  (*	      2.  General denominator section to a symbolic power	*)
  (*		  Step is different because we don't know what k is.	*)
  myinvz[ c_. / (1 + b_. z_^-1)^k_., z_, n_, rm_, rp_, rest___ ] :>
	leftOrRightSided[b, c Binomial[n+k-1,k-1] (-b)^n Step[n+k-1],
			 n, rm, rp] /;
	FreeQ[{b,c,k}, z],

  (*	      3.  General numerator section to a symbolic power		*)
  myinvz[ c_. (a_ + b_. z_^-1)^k_, z_, n_, rm_, rp_, s_, op_, rest___ ] :>
	Block [ {},
		Assuming[ Positive[k], dialogueAllFlag ];
		c leftOrRightSided[ -b/a, Binomial[k,n] b^n a^(k-n) Step[n],
				    n, rm, rp ] ] /;
	FreeQ[{a,b,c,k}, z] && Implies[NumberQ[k], N[k > 0]],

  (*	      4.  Reduce expressions in z^k to expressions in z.	*)
  (*		  This applies to expressions other than polys in z.	*)
  myinvz[ f_, z_, n_, rest___ ] :>
	upsampledZtransform[ f, z, n, rest ] /;
	ZUpsampledQ[f, z],


  (*	  E.  Factor rational polynomials and make the denominator	*)
  (*		only contain negative powers of z.			*)

  (*	      1.  Factor rational polynomials when denominators have	*)
  (*		  positive powers of z (ensured by PolynomialQ test).	*)
  myinvz[ p_, z_, n_, rm_, rp_, s_, rest___ ] :>
	myinvz [ ( Expand[Numerator[p]] ) / ( Factor[Denominator[p]] ),
		 z, n, rm, rp,
		 SetStateField[s, polyfactorfield, False], rest ] /;
	GetStateField[s, polyfactorfield] && PolynomialQ[Denominator[p], z] &&
	  RationalFunctionQ[p, z],

  (*	      2.  Make factored denominators have neg. powers of z	*)
  (*		  a.  For numeric k, numeric j				*)
  myinvz[ ((a_ + b_. z_^j_.)^k_) x_., z_, n_, rest___ ] :>
	b^k myinvz[ Distribute[Numerator[x] z^(j k)] ( 1 + a z^(-j) / b )^k /
		    Denominator[x], z, n, rest ] /;
	FreeQ[{a,b}, z] && IntegerQ[k] && (k < 0) && IntegerQ[j] && (j > 0),

  (*		  b.  For symbolic k, numeric j				*)
  myinvz[ x_. / ((a_ + b_. z_^j_.)^k_), z_, n_, rest___ ] :>
	b^-k myinvz[ Distribute[Numerator[x] z^-(j k)] / (1 + a z^(-j) / b)^k /
		     Denominator[x], z, n, rest ] /;
	FreeQ[{a,b,k}, z] && IntegerQ[j] && (j > 0),

  (*	      3.  Normalize factored denominators.			*)
  (*		  a.  For numeric k, numeric j				*)
  myinvz[ ((a_ + b_. z_^j_)^k_) x_., z_, n_, rest___ ] :>
	a^k myinvz[ x ( 1 + b z^j / a ) ^ k, z, n, rest ] /;
	FreeQ[{a,b}, z] && IntegerQ[k] && (k < 0) && IntegerQ[j] &&
	  (j < 0) && ! SameQ[a, 1],

  (*		  b.  For symbolic k, numeric j				*)
  myinvz[ x_. / ((a_ + b_. z_^j_)^k_), z_, n_, rest___ ] :>
	a^-k myinvz[ x / ( 1 + b z^j / a ) ^ k, z, n, rest ] /;
	FreeQ[{a,b,k}, z] && IntegerQ[j] && (j < 0) && ! SameQ[a, 1],




  (*  II. N O N - R A T I O N A L     T R A N S F O R M S  *)


  (*      A.  Handle Exp[affine] forms; covers Sinh[affine]	*)
  (*          Cosh[affine] forms converges for any ROC		*) 
  myinvz[ c_^(b_. + a_. z_^-1), z_, n_, rm_, rp_, rest___ ] :>
	leftOrRightSided[ a Log[c], c^b a^n (Log[c])^n Step[n] / n!,
			  n, rm, rp ]  /;
	FreeQ[{a,b,c}, z],


  (*	  B.  Handle Log[affine] forms [O&S, 55]		*)
  myinvz[ Log[a_ + b_. z_^-1], z_, n_, rm_, rp_, rest___ ] :>
	Log[a] Impulse[n] +
	  leftOrRightSided[-b/a, -(-b/a)^n Step[n-1] / n, n, rm, rp] /;
	FreeQ[{a,b}, z],


  (*	  C.  Family of causal IIR filters			*)
  myinvz[ ( 1 + z_^-1 )^r_, z_, n_, rm_, rp_, s_, op_, rest___ ] :>
	Block [ {},
		Assuming[ Positive[r + 1], dialogueAllFlag ];
		Gamma[1 + r] Step[n] / ( Gamma[1 + n] Gamma[1 + r - n] ) ] /;
	FreeQ[r,z] && Implies[ NumberQ[N[r]], N[r > -1] ],


  (*	  D.  Other unusual forms.				*)
  myinvz[ ArcTan[z_^-1], z_, n_, rest___ ] :>
	Sin[n Pi / 2] Step[n - 1] / n,

  myinvz[ Sqrt[z_] ArcTan[Sqrt[z_^-1]], z_, n_, rest___ ] :>
	Step[n] / ( 2 n + 1 ),

  myinvz[ Sqrt[z_] ArcTan[z_^-1], z_, n_, rest___ ] :>
	Step[n] / ( 2 n + 1 ),

  myinvz[ Cosh[ Sqrt[z_^-1] ], z_, n_, rest___ ] :>
	Step[n] / ( 2 n ) !,

  myinvz[ Sinh[ Sqrt[z_^-1] ], z_, n_, rest___ ] :>
	Step[n] / ( 2 n + 1 ) !,

  myinvz[ Sqrt[ 1 / ( 1 + a_. z_^-1 ) ], z_, n_, rest___ ] :>
	(-a)^n (2 n)! Step[n] / ((2^n n!) ^ 2)  /;
	FreeQ[a, z],




  (* III. Z - T R A N S F O R M     P R O P E R T I E S			*)

  (*  Omitted these properties:	1.  multiplication by n^m		*)
  (*				2.  multiplication by Cos[bn]		*)
  (*				3.  multiplication by Sin[bn]		*)
  (*				4.  multiplication by x y		*)
  (*									*)
  (*  This section is similar to the z-transform Properties section	*)
  (*  of the forward z-transform rules.					*)


  (*	  A.  Downsampled sequence.					*)
  (*          Had to put it before mult. by exponential sequence	*)
  myinvz[ Summation[k_, 0, Mminus1_, 1][f_], z_, n_, rm_, rp_, rest___ ] :>
	Block [	{m, newf, newrm, newrp},
		m = Mminus1 + 1;
		newf = downsampledRewrite[f, k, m, z];
		newrm = rootROC[rm, 1/m];
		newrp = rootROC[rp, 1/m];
		m Downsample[m, n][myinvz[newf, z, n, newrm, newrp, rest]] ] /;
	downsampledQ[f, k, Mminus1 + 1, z],


  (*	  B.  Multiplication by exponential sequence (exp)^n		*)
  (*	      Check this possibility out before partial fractions	*)
  (*	      Moved from rational section because it interfered		*)
  (*	      with the Log lookup.  This rule should not interfere	*)
  (*          with the downsampling rule which requires a symbolic	*)
  (*	      summation to be present.					*)
  myinvz[ f_, z_, n_, rm_, rp_, s_, rest___ ] :>
	Block [	{exp, fun, newf, newrm, newrp, state},
		exp = 1 / ScalingFactor[f, z];
		newf = f /. z -> z exp;
		newrm = Abs[rm / exp];
		newrp = If [ InfinityQ[rp], Infinity, Abs[rp / exp] ];
		state = SetStateField[s, expcheckfield, False];
		(exp)^n myinvz[newf, z, n, newrm, newrp, state, rest] ] /;
	GetStateField[s, expcheckfield] && FreeQ[f, Summation] &&
	multByExponentialQ[f, z],


  (*	  C.  Homogeneity --  pick off constants			*)
  myinvz[ c_ x_, z_, n_, rest___ ] :>
	c myinvz[ x, z, n, rest ] /;
	FreeQ[c, z],


  (*	  D.  Delay or shift						*)
  myinvz[ x_ z_^m_., z_, n_, rest___ ] :>
	shift[ myinvz[x, z, n, rest], n, -m ] /;
	FreeQ[m, z] && Implies[ NumberQ[m], RationalQ[m] ],


  (*	  E.  Simple Additivity						*)
  myinvz[ x_ + y_, z_, n_, rm_, rp_, s_, rest___ ] :>
	myinvz[x, z, n, rm, rp, resetstate[s], rest] +
	myinvz[y, z, n, rm, rp, resetstate[s], rest],


  (*      F.  Partial fractions decomposition				*)
  (*	      partial fractions decomposition for denominators that	*)
  (*	      are polynomials in of z; once partial fractions de-	*)
  (*	      composition has been applied to an expression, it is	*)
  (*	      never applied to any of the resulting sub-expressions.	*)
  (*	      This is controlled by the partial fractions semaphore/	*)
  (*	      flag in the local state; at this point, the denominator	*)
  (*	      should be a function of z^-1, but we will do partial	*)
  (*	      fractions on a denominator that is a "two-sided" poly-	*)
  (*	      nomial (one in both positive and negative	powers of z).	*)

  (*	      1.  Normalize the denominator				*)
  myinvz[ p_, z_, n_, rm_, rp_, s_, rest___ ] :>
	myinvz[ normalizeDenominator[p, z], z, n, rm, rp,
		SetStateField[s, normdenomfield, False], rest ] /;
	GetStateField[s, normdenomfield] && partialFractionsQ[p, z],

  (*	      2.  Use Mathematica's Apart primitive			*)
  myinvz[ p_, z_, n_, rm_, rp_, s_, rest___ ] :>
	Block [ {failure, partfrac, pnorm, pfnorm, state},
		partfrac = KeepNormalized[p, Apart, z, z];
		state = SetStateField[s, partialfractionsfield, False];
		failure = partialFractionsQ[partfrac, z] &&
		          ( SameQ[partfrac, p] ||
			    SameQ[normalizeDenominator[p, z],
				  normalizeDenominator[partfrac, z]] );
		If [ failure,
		     partfrac = p,
		     partialFractionsDialogue[p, partfrac, op, Apart] ];
		myinvz [ partfrac, z, n, rm, rp, state, rest ] ] /;
	GetStateField[s, partialfractionsfield] && partialFractionsQ[p, z],

  (*	       3.  Do the partial fractions that Apart cannot do	*)
  myinvz[ p_, z_, n_, rm_, rp_, s_, op_, rest___ ] :>
	Block [ {filter, partfrac, state},
		filter = If [ SameQ[Replace[Apart, op], All], N, Identity ];
		partfrac = KeepNormalized[p, MyApart, z, z, filter];
		state = SetStateField[s, partialfractionsfield, False];
		state = SetStateField[state, partialfractionsfield2, False];
		If [ ! SameQ[partfrac, p],
		     partialFractionsDialogue[p, partfrac, op, MyApart] ];
		myinvz [ partfrac, z, n, rm, rp, state, op, rest ] ] /;
	GetStateField[s, partialfractionsfield2] && partialFractionsQ[p, z],


  (*	  G.  Additivity of Numerator					*)
  myinvz[ (x_+y_)/c_, z_, n_, rm_, rp_, s_, rest___ ] :>
	myinvz[x/c, z, n, rm, rp, resetstate[s], rest] +
	myinvz[y/c, z, n, rm, rp, resetstate[s], rest],


  (*	  H.  Logarithm properties 					*)
  myinvz[ Log[a_. / b_], z_, n_, rest___ ] :>
	myinvz[Log[a] - Log[b], z, n, rest],

  myinvz[ Log[a_ b_], z_, n_, rest___ ] :>
	myinvz[Log[a] + Log[b], z, n, rest],


  (*      I.  Derivative						*)
  myinvz[ Derivative[m_][f_][z_], z_, n_, rest___ ] :>
	shift[ ZPolynomial[m, n] myinvz[f[z], z, n, rest], n, -m ] /;
	FreeQ[m, z],




  (* IV.  S I G N A L     P R O C E S S I N G     S T R U C T U R E S	*)


  (*        -.  An operator independent of n---  take z-transform of f    *)
  myinvz[ operator_[params__][f_], z_, n_, rest___ ] :>
	operator[params][ myinvz[f, z, n, rest] ] /;
	FreeQ[{params}, z],


  (*	  A.  Difference equations					*)
  myinvz[ x_ (1 - z_^-1)^m_., z_, n_, rm_, rp_, s_, op_, rest___ ] :>
	Block [ {},
		Assuming[ m < 0, dialogueAllFlag ];
		If [ dialogueAllFlag, Print[ "where ", m, " is an integer" ] ];
		Difference[m,n][ myinvz[x, z, n, rm, rp, s, op, rest] ] ] /;
	FreeQ[m,n] && Implies[NumberQ[m], IntegerQ[m] && m > 0],


  (*	  B.  Summation							*)
  myinvz[ x_ / ( 1 - z_^-1 ), z_, n_, rest___ ] :>
	Block [	{context, index},
		context = $Context;
		$Context = "Global`";
		index = Unique["index"];
		$Context = context;
		Summation[index,0,n,1][ myinvz[x, z, index, rest] ] ],


  (*	  C.  Periodic sequence with period k.				*)
  myinvz[ x_ / ( 1 - z_^k_), z_, n_, rm_, rp_, s_, op_, rest___ ] :>
	Block [ {},
		Assuming[ k < 0, dialogueAllFlag ];
		If [ dialogueAllFlag, Print[ "where ", k, " is an integer" ] ];
		Periodic[-k,n][ myinvz[x, z, n, rm, rp, s, op, rest] ] ]  /;
	FreeQ[k,n] && Implies[ NumberQ[k], IntegerQ[k] && ( k < 0 ) ],




  (*  V.  Z - T R A N S F O R M    S T R A T E G I E S			*)

  (*        If the following three rules are considered, then the	*)
  (*  the current expression could not be inverted by table lookup	*)
  (*  and applying properties.  So, the following strategy is used:	*)

  (*      A.  Factor terms inside of a log or exponential function.	*)
  (*      B.  Normalize the denominator.				*)
  (*      C.  Expand every term in the z-transform to its simplest	*)
  (*	      components and try again.					*)
  (*      D.  Substitute 1/z for z, take the inverse z-transform, and	*)
  (*	      substitute -n for n in the result.			*)
  (*	  E.  Complex cepstrum.						*)
  (*	  F.  If the current expression is a rational polynomial,	*)
  (*	      then return an ARMA signal processing expression since	*)
  (*	      it could not be inverted using partial fractions.		*)
  (*      G.  Perform a Laurent series expansion using the first	*)
  (*	      Terms terms (Terms is an option) and try again.		*)


  (*      A.  Factor terms inside of a logarithmic or exponential expr.	*)
  myinvz[ Log[x_], z_, n_, rm_, rp_, s_, rest___ ] :>
	myinvz [ Log[Factor[x]], z, n,
		 rm, rp, SetStateField[s, logformfield, False], rest ] /;
	GetStateField[s, logformfield] && PolynomialQ[x, z],

  myinvz[ x_. b_^(exp_), z_, n_, rm_, rp_, s_, rest___ ] :>
	myinvz [ x b^(Factor[exp]), z, n,
		 rm, rp, SetStateField[s, expformfield, False], rest ] /;
	GetStateField[s, expformfield] && FreeQ[b, z],


  (*      B.  Expand definitions of all expressions to simplest form	*)
  myinvz[ x_, z_, n_, rm_, rp_, s_, rest___ ] :>
	myinvz [ ExpandAll[x], z, n, rm, rp,
		 SetStateField[s, expandallfield, False], rest ] /;
	GetStateField[s, expandallfield],


  (*      C.  Substitute z^-1 for z, take the inverse z-transform,	*)
  (*	      and substitute -n for n in the result.			*)
  myinvz[ x_, z_, n_, rm_, rp_, s_, rest___ ] :>
	Block [	{newrm, newrp, newx, seq},
		newx = x /. z -> z^-1;
		newrp = If [ ZeroQ[rm], Infinity, 1/rm ];
		newrm = If [ InfinityQ[rp], 0, 1/rp ];

		seq = MyInvZTransform [ newx, z, n, newrm, newrp,
					anticausalInvZstate[], rest ];
		If [ InvalidInvZTransformQ[seq],
		     myinvz [ x, z, n, rm, rp,
			      SetStateField[s, anticausalfield, False], rest ],
		     seq /. n -> -n ] ] /;
	GetStateField[s, anticausalfield],


  (*	  D.  Complex cepstrum [O&S, 498]				*)
  myinvz[ Log[x_], z_, n_, rm_, rp_, s_, rest___ ] :>
	-1/n myinvz [ z D[x, z] / x, z, n, rm, rp,
		      SetStateField[s, cepstrumfield, False], rest ] /;
	GetStateField[s, cepstrumfield] &&
	  FreeQ[ D[x, z], Derivative[1][f_][z] ],


  (*	  E.  Invert rational polynomial as a cascade of an		*)
  (*	      FIR filter and an IIR filter.				*)
  myinvz[ x_, z_, n_, rm_, rp_, s_, rest___ ] :>
	Block [	{firlist, iirlist, maxexp, zdenom, znumer},
		firlist = CoefficientList[znumer /. z -> z^-1, z];
		iirlist = CoefficientList[zdenom /. z -> z^-1, z];
		If [ SameQ[znumer, 1],
		     Shift[ maxexp, n ]
			  [ IIR[n, iirlist] ],
		     Shift[ maxexp, n ]
			  [ FIR[n, firlist]
			    [ IIR[n, iirlist] ] ] ] ] /;
	RationalFunctionQ[x, z] && RationalPolynomialQ[x, z]

}


(*  E N D     P A C K A G E  *)

End[]
EndPackage[]

If [ TrueQ[ $VersionNumber >= 2.0 ],
     On[ General::spell ];
     On[ General::spell1 ] ];


(*  H E L P     I N F O R M A T I O N  *)

Combine[ SPfunctions, {InvZTransform} ]
Protect[ InvZTransform ]


(*  E N D I N G     M E S S A G E  *)

Print["The inverse z-transform rules have been loaded."]
Null
